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r^ ' Abstract 

(N ! 

We consider the problem of computing the distance between two piecewise-hnear bivariate 
^ ' functions / and g defined over a common domain M. We focus on the distance induced by 

the i2-norm, that is ||/ — g\\2 ~ y IlAiif " d)^ ■ ^^ f ^^ defined by hnear interpolation over 

a triangulation of M with n triangles, while g is defined over another such triangulation, the 

obvious naive algorithm requires 8(ri,^) arithmetic operations to compute this distance. We 

. show that it is possible to compute it in 0{nlog n) arithmetic operations, by reducing the 

\i^ ' problem to multi-point evaluation of a certain type of polynomials. 

^^ . We also present an application to terrain matching. 
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I Introduction and problem statement 

In this paper we use a novel combination of tools from computational geometry and computer 
algebra to speed up a computation of || • ||2-norm distance between two bivariate piecewise-linear 
functions. Algebraic tools have already been used in other recent work in computational geometry, 
to seemingly defy "obvious" lower bounds: for example, Ajwani, Ray, Seidel, and Tiwary [2] use 
algebraic tools to compute the centroid of all vertices in an arrangement of n lines in the plane, 
without explicitly computing the vertices. 

We feel that working on computational geometry problems using a combination of traditional 
and algebraic tools expands the repertoire of questions that can be approached and answered 
satisfactorily. Employing such combinations of methods expands the horizon of solvable problems. 
Indeed, in a significant recent development, several breakthrough results have been obtained by 
applying algebraic methods to problems of combinatorial geometry: Guth and Katz's recent work 
on the joints problem |8j and on the Erdos distinct distance problem [9] has triggered an avalanche 
of activity; see, for example, [5l[6| [T0l - [l4] . It appears that the use of algebraic tools in geometry 
(both combinatorial and computational) allows one to approach problems inaccessible by more 
traditional methods. The present work is just one step in that direction. 

Background and previous ■work In [I], Aronov et al. considered a quite common object in 
computational geometry and geographical information systems (GIS): that of a "terrain." A terrain 
is (the graph of) a bivariate function over some planar domain, say, a rectangle or a square. It 
is often used to model geographic terrains, e.g., elevation in a mountainous locale, but can also 
be applied to storing any two-dimensional data sets, such as precipitation or snow cover data. 
A common (though by no means the only) way of interpolating and representing discrete two- 
dimensional data is a triangulated irregular network {TIN): the values of a bivariate function are 
given at discrete points in, say, the unit square. The square is triangulated, using data points as 
vertices. The function is then linearly interpolated over each triangle. This produces a piecewise- 
linear approximation of the "real" (and unknown) function. 

The problem raised in [1] was that of comparing two terrains over the same domain, say, the unit 
square, but given over two unrelated triangulations. One could imagine comparing the outcome of 
two different ways of measuring the same data, or finding correlation between, say, the elevation 
and the snow cover over the same geographic region. The focus of that work was on identifying 
linear dependence between the two functions or terrains. Three natural distance measures between 
the two functions were considered and several algorithms presented for computing such a distance 
and optimizing it, subject to vertical translation and scaling. The only observation made for the 

II • II2 norm (see the definitions below) in [1] is that, if the two terrains share a triangulation, both 
the distance computation and the optimization problem can be solved easily in linear time, while 
it appears that in general quadratic time seems to be needed to deal with the case of arbitrarily 
overlapping triangulations. The substance of the current work is disproving this assertion and 
describing a near-linear-time algorithm for both problems. 

Problem statement and results Given bivariate functions f,g: M — )• M, one can naturally 
define a distance between them as 

\\f-9h = {jj {fix, y) - g{x, y)fdxdy) ^'\ 

Expanding the expression under the integral, we obtain //(/ — 5)^ = Jf f^ — 2 Jf fg + J g'^. If the 
two functions are piecewise linear, defined over different triangulations of M, only the middle term 



presents a problem for efficient computation. Thus, in the bulk of the paper we will focus on the 
computation of JJ fg, showing the following: 

Theorem 1.1. Given piecewise-linear functions f and g defined over different triangulations of the 
same domain M, with n triangles each, fjj^,jf{x,y)g{x,y)dxdy can be computed using (9(?ilog n) 
arithmetic operations. 



Armed with this result, as already mentioned, we can quickly compute ||/ — g\ 
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Theorem 1.2. Given piecewise-linear functions f and g defined over different triangulations of 
the same domain M, with n triangles each, \\f — g\\2 can be computed using 0(nlog n) arithmetic 
operations. 

Naively, the integral in Theorem 11.11 can be expressed as a sum of integrals over each cell 
appearing in the overlay of the two triangulations of / and g. Unfortunately, this overlay has a 
quadratic number of cells, in the worst case. The main idea of our algorithm is to reduce the 
computation of the integral to double sums of algebraic functions over grids, which allows us tu 
use fast multi-point evaluation algorithms. 

The paper is organized as follows. In section [21 we will show how the integral of a function 
over a convex polygon can be expressed as a sum of elementary algebraic functions over its vertices. 
Applying the process to a convex decomposition of a region M expresses the value of JJ^j fg as a 
summation over the vertices of the decomposition. Then, in section [3l we show how the integral 
over the overlay of the two triangulations can be reduced to a sum of elementary functions over 
pairs of edges, plus some additional terms computable in linear time. Li section [U we use the 
bipartite clique decomposition to arrange the pairs of edges in complete grids of fairly rigid form. 
Finally, in section [U we use a fast multi-point evaluation algorithm to compute the sums over each 
grid, completing the description of our method. 

2 How to integrate over a convex subdivision 

Consider a bivariate function h defined over a convex polygon C in the plane. To simplify our 
presentation and without loss of generality, we assume that all vertices of C lie to the right of the 
y-axis. For a vertex p of C we refer to the lines supporting the edges of C incident to p as L{p, C) 
(the one with higher slope) and U{p, C) (the one with lower slope); to the left of p, L(p, C) is below 
U{p,C). Let 

L{P, C): y = yi{p, C) + si{p, C)x, and 
U{p, C):y = yu{p, C) + s„(p, C)x. 

We omit the explicit dependence on C and/or p whenever it causes no confusion. Define 
6{p,C)-- 



— 1 if C is above both L{p, C) and C/(p, C), or below both of them, 
+1 if C is below U{p, C) and above U{p, C), or vice versa. 



Finally, put 



/■^Pj ryu{p,C)+su{p,C)x 
T{p,C,h) := 6{p,C) / h{x,y)dxdy. 

Jx=0 Jy=yi{p,C)+si{p,C)x 



In words, T is the signed integral of h over the triangle Tp delimited by the y-axis, L{p,C), and 
U{p,C). 




Figure 1: An illustration of the proof of Lemma 12. li Tr is dotted; T^ is dashed; Ax is lightly 
shaded. 

With the above notation, we express the integral of h over C in a convenient way as a sum of 
terms associated with its vertices: 

Lemma 2.1. Let C be a convex polygon with vertices pi, ■■■,Pk, and h a bivariate function. Then 
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h{x, y)dxdy = ^ T{pj,C, h). 



(1) 



Proof. Partition the vertices of C into four subsets Vl, Vr-, Vt, and Vb as follows. Vl := {pl} 
(resp., Vr := {pr}) consists of the unique leftmost (resp., rightmost) vertex of C. Vr := {ti, • • •} 
(resp., Vb := {bi, . . .}) is the sequence of vertices of C from p^ to pi in the counterclockwise (resp., 
clockwise) direction; refer to Fig. [T] 

To each vertex p, we associate the triangle Tp as defined above. Put Tl := Tp^ and Tr := Tp^^. 
Since C is convex, so for all i, line tjtj+i (resp., bibi+i) is below the line tj-itj (resp., above the 
line bi-ibi) left of C. Thus, the triangles Tj, t € Vr, (resp., Tb, for b € Vr) do not overlap. Let 
At = UtGV ^p ^11"^ ^B = \JbeV '^b- By construction we have 

AtVJAb = {TlVJTr)\C and ^t n ^b = T^ n T^. 

Since C C TlL) T^i, the first equality, written in terms of characteristic functions, gives 

Iat + l^s — IathAs = Iti, + Ith - ^TlCiTr — lc> 
which can be simplified, using the second equality, to yield, as promised 



IAt + IAs = iTr + ITp - Ic- 
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We now consider a convex subdivision C of some bounded region M in the plane, with each cell C 
associated with its own bivariate function he, thereby defining a function h over all of M (as it does 
not affect the value of the integral, the functions he need not agree along the common boundaries 
of adjacent cells). By summing eq. ([T|) over the cells C of C, we can compute ffj^j h{x, y)dxdy: 



Corollary 2.2. 



h{x,y)dxdy= ^ j l^hc{x,y)dxdy = ^^ T{p,C,hc). (2) 

"" cell c ec 

C adjacent to p 



^ cell C &C'' ''^ cell C €C 



One of the advantages of the formulation in eq. ([2]) for our application is that an individual 
integral under the sum can be expressed as a rational function when he is a bivariate polynomial. 

Lemma 2.3. Let p be an intersection point of y = yi + six and y = yu + s^x as above. Then 

_ Vu -yi 

Xp — 

Su- Si 

and 

x=0 Jy=yi+six W« ^i) 

where Pij is a polynomial of total degree i + 2j + 2. 

Proof. Let Qij{u,v,x) be the only polynomial such that —q^ = x*(u + vxy and Qij{u,v,0) = 
for all u,v &W. In particular, Qij has total degree i + 2j + 1, its degree in x is i + j + 1, and the 
coefficient of x'"'"-'"'""^ in Qij is j^rj^v^ . Now 

/ x'yHxdy = / —-—{x\yu + s«x)-''+^ - x\yi + sixy^^)dx 

x=Q J y=yi+six J x=0 J ' ^ 

=——r{Qi,j+i{yu,Su,Xp) - Qij+i{yi,si,Xp)). 

Therefore the value of the integral can be expressed as a polynomial in yi,yu,si,Su,Xp of total 
degree i + 2j ' + 3 of the form 



1 j+i+i 

—jist^ - si^^)4^^^^ + YI 1k{yu,Su,yusi)xl. 



fc=0 

After substituting — ^""^' for Xp and bringing to a common denominator, we conclude that the 
expression can be rewritten in the form 

{su- si)P{yi,yu,si,Su) 

{su - siy+^+^ 

as claimed. D 

3 How to integrate over an overlay 

In our problem, we are interested in computing the integral JJ f{x, y)g{x, y)dxdy, where / is defined 
over a triangulation Tj of M C M?, with a separate linear function fA{x,y) determining / over 
each triangle A € Tj; g is defined similarly over a different triangulation Tg of M. The product 
h{x, y) := f{x, y)g{x, y) is thus naturally defined over the convex decomposition C of M that is the 
overlay of Tj and T^. By Corollarv 12.21 it is sufficient to evaluate a sum over all vertices of C The 



vertices of C come in two flavors: the original vertices of Tj and of Tg, and the intersections of 
edges of Tj and T^. Therefore, 

f{x,y)g{x,y)dxdy = ^^ T{p,C,hc) 

p vertex of C 
C adjacent to p 

Y, Y. np,c,hc)+ Y E np,c,hc). 

p vertex C adjacent p=eine2, C adjacent 

(ei,e2)GT/xT9 top 



'Ey involves 0{n) integrals. We preprocess each of Tj and T^ for logarithmic-time point location 
queries, in 0{nlogn) time (see, for example, [3]). For each vertex p of Tj, we locate the triangle 
Ag G Tg containing it. Then, for each triangle Aj £ Tj incident to p, we can compute T{p, Ag n 
Ah, fAfdAg) in constant time. The treatment of vertices of T^ is symmetric. We spend 0{nlogn) 
total for point location. The remaining work is proportional to the sum of vertex degrees in both 
triangulations, which is 0{n). Hence S^ can be computed in 0{nlogn) arithmetic operations. We 
devote the rest of the discussion to computing Eg efficiently. 

4 Bipartite clique decomposition 

Let Ef be the set of edges of T/- and Eg the set of edges of Tg. In [4], it is shown that it is possible 
to compute a family T = {{Ri, Bi), . . . , {Ru, Bu)} where Rk C Ef and Bk C Eg, such that 

(i) every segment in Rk intersects every segment in Bk] 

(ii) every segment of R^ has lower slope than every segment of B^, or vice versa; 
(iii) for every intersecting pair (ei, 62) £ Ej x Eg there is exactly one k such that ei £ Rk, €2 £ Bk] 

no such k exists for a non-intersecting pair (ei, 62) £ Ef x Eg] 
(iv) Eki\Rk\ + \Bk\) = Oinlog^n). 

This family can be computed in 0{nlog n) time. 

5 Multi-point evaluation 

In this section, we explain how to efficiently compute 

Y E r(p,c,M=EE E np,c,hc), (3) 

p=eine2, C adjacent eig/? 626-6 C adjacent 

{ei,e2)€RxB top to p = ei n 62 

where {R, B) := {Rk,Bk) is a pair of sets of triangulation edges produced by the bipartite clique 
decomposition. In particular, all segments of R intersect all segments of B and, moreover, without 
loss of generality, the slopes of segments of R are greater than those of segments of B. 

5.1 Reduction to sums of rational functions 

Each triangle of Tf and Tg is associated with a bivariate linear function. If e is an edge of Tj U T^, 
let fu{e) be the linear function associated to the upper triangle (in the triangulation to which e 
belongs) adjacent to e; fi{e) is the corresponding function for the lower triangle; we can define the 
function to be identically zero for regions outside M, but it will never be used by the algorithm. 



A vertex p = ei n 62, with ei £ R and 62 € B, lies on the boundary of four cehs of C. We focus 
on the cell Cieft = Cieft{p) lying above ei and below 62, for which p is the rightmost point. Thus 
^Ci^ft = fu{ei)fiie2) and 6{p,Cieft) = +1- Suppose /«(ei): (x,y) ^ a(ei) + 6(ei)x + c(ei)y and 
fi{e2) : (2;, y) i-> 0(62) + &(e2)a; + c{e2)y- We compute the contribution to the sum ([3]) of such cells, 
over all choices of p. (The remaining three types of cells adjacent to p are treated by an entirely 
symmetric argument.) Given an edge e of Tj U T^, let y = y{e) + s{e)x be the equation of the line 
supporting it, so we can write 



T{p, Cleft, fu{ei)fu{e2)) = 

i-Xp ry{e2)+s(e2)x 
Jx=0 Jy=y{ei)+s{ei)x 



a(ei)a(e2) + (a(ei)6(e2) + a{e2)b{ei))x + 
(a(ei)c(e2) + a(e2)c(ei))y + b{ei)b{e2)x'^ + | dxdy. 
{b{ei)c{e2) + b{e2)c{ei))xy + c{ei)c{e2)y'^ 



The above integral can be expressed as the sum of nine integrals of a function of the form 
v{ei)'w{e2)x^y^ , where v and w are some functions that assign a real number to each edge. Gather- 
ing all the terms and recalling that p = ei n 62 , eq. ([3]) can be rewritten as the sum of 36 expressions 
of the form 



eieRe2GB 



a;eine2 rvie2)+s{e2)x 

/ v{ei)w{e2)x^y^ 

x=0 J y=y{ei)+s{ei)x 



v^ v^ v{ei)w{e2)Pi,j{y{ei),y{e2), s{e i), 3(62)) 
^ ^ (s(e2) - Kei))^+i+i 



v^ v\(^\)'w\(^2)n,jKyK(i\),yKfi2),sKex),s\e2)) , . 

eii=Re2eB 



by Lemma 12.31 

5.2 Fast multi-point evaluation 

Now we will use multi-point evaluation to speed up the computation of eq. (JH). To accomplish this, 
we will replace the values associated to the edges of R by symbolic variables, while using the actual 
numerical values for those for the edges of B. Then we will compute the corresponding symbolic 
rational function using a divide- and-conquer strategy (Lemma lS.ip . Finally, we will use multi-point 
evaluation on the resulting polynomials (Lemma 15.21) . 






eS-B 

can he expressed in the form 



Lemma 5.1. Let u and v be two functions from the edges of B to M and suppose \B\ < n. Then 

u{e) 
{X-v{e)Y 

D[xy 

where N{X) and D{X) are polynomials of degree at most (n — l)d and nd respectively; their 
coefficients can he computed explicitly in 0{J^{nd)\ogn) arithmetic operations, where Ad{q) = 
0{q\ogq) is the cost of multiplication of two univariate polynomials of degree at most q. 

Proof. For simplicity of presentation and without loss of generality, assume that \B\ is a power of 
two. We bring eq. (I5|) to a common denominator by combining the fractions in pairs, reducing their 

6 



number to \B\/2, and repeating the process log \B\ = ©(logn) times. The bounds on the degree of 
the final numerator and denominator are immediate from examining the original fractions. 
We now explain how to bring 

Ni{X) N2{X) 

with Ni,N2,Di,D2 of degree at most kd, to a common denominator in time 3A4{kd) + 0{kd). 
Indeed, the above fraction is equal to 

Ni{X)D2iX) + N2iX)Di{X) 
Di{X)D2{X) 

so it can be computed by three calls to fast polynomial multiplication plus a linear number of 
additional operations, as claimed. 

This completes the proof of the lemma, as the cost of one round of combining fractions with 
denominators and numerators of degree at most kd is n/k ■ {2>M.{kd) + 0{kd)) = 0{A4{nd)), since 
M is superlinear. D 

The second lemma handles summing the values of a special kind of polynomials. 

Lemma 5.2. Let P{Xq,Xi, . . . ,Xr) be a polynomial of degree n in Xq and d in Xi, . . . ,Xr. Let 
E be a set of at most n points ofW^^. The values of P at the points of E can be simultaneously 
computed in 0(( ^'') A^ (n) log n) time. 

Proof. P can be expanded with respect to the variables Xi , . . . , X^ , and has at most ( ^^) mono- 
mials. Each coefficient is a univariate polynomial in Xq of degree at most n. Then, using standard 
multi-point evaluation algorithm for univariate polynomials [3 ch. 10], we can compute simulta- 
neously the values of these univariate coefficients of P in 0(( ^)A4{n) logn) at every point of E. 
Finally, we compute the values of each of the monomials of P in time ©(( ^^^n). Combining all 
these values also costs 0(( )n) arithmetic operations, concluding the proof. D 

Now we are ready to efficiently evaluate eq. ([H). Let F{X,Y,V) be the polynomial 

After expanding the numerator of this fraction, we note that F has the form 

{s{e2)-X) 



F{x,Y,v) = E ( E ''"''";!;i"!^!'^'l?l'i'^''^^ ) x'^-y'^-v. 



0<dx<i+2j+2 \e2GS 
0<dY<i+2j+2 

In our case, i + j + 1 < 3, and, using Lemma l5.H we can compute each coefficient of X^Y^V 
and express F in the following form, in 0{nlog n) time: 

^ ' '^ DiX) ' 

with N a polynomial of degree n{i + j + 1) + 1 in X, i+ j + 2 in Y, and D a univariate polynomial 
of degree n{i + j + 1). 



Finally, eq. (j3|) can be rewritten as 

v^ v^ viei)w{e2)Pi,j{y{ei), y{e2), s{ei), s{e2)) ^ >p N{s{ei),y{ei))v{ei) 
^n^B {s{e2) - sie,)y-^-^ ^^^ I.(.(e,)) " 

As \R\ < n, using Lemma 15.21 we can compute simultaneously all the terms under the sum, and 
add them together in ©(nlog n) arithmetic operations. 

5.3 Putting it together 

To summarize, we can evaluate eq. (jl]) in 0((|i?fc| + |i?jt|)log {\Bk\ + |i?A;|)) operations, for each 
pair {Rk, Bk). We also saw in section d] that X^^d-Rfel + l^fel) = 0{n\og^ n), where n is the number 
of triangles appearing in our triangulations. Therefore 

J2i\Rk\ + \Bk\)log\\Rk\ + \Bf,\) < Y.{\Rk\ + \Bk\)log^n = ©(nlog^n), 
k k 

which allows us to conclude that we can compute Eg in ©(nlog n) time. This completes the proof 
of Theorem II. 1[ 

6 Discussion 

In [1], the following optimization problem was considered: given two functions / and g and a 
distance measure ||/ — ffH between them (the paper discusses || • ||i, || • ||2, and || • ||oo, but we only 
consider || -112 here), find the values of real parameters s and t that minimize ||/ — (sg+i)!!. If/and(7 
are interpreted as geometric "terrains," we are looking for the scaling and translation of the vertical 
coordinate of the terrain g to best match terrain / [l]. Since ||/ — {sg + t)\\\ = //(/ — {sg + 1))^ is 
a degree- two polynomial in s and t with coefficients easily expressible in terms of JJ /, Jjg, fff'^, 
ffg^, Jf fg, and fj 1, being able to compute Jf fg efficiently immediately yields 

Theorem 6.1. Given piecewise-linear functions f and g defined over different triangulations of 
the same domain M, with n triangles each, the values s and t minimizing \\f — {sg + t)||2 can be 
computed using 0(nlog ?i) arithmetic operations. 

Theorems 11.11 II. 2( and 16.11 extend to piecewise-polynomial functions of constant maximum 
degree with essentially no modifications. What other classes of functions can be handles using 
similar methods? 
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